学习自 https://confluence.cornell.edu/pages/viewpage.action?pageId=333371302

以ANSYS 17.0为例。

该算例分为两个部分，第一部分将采用运动参考系（Moving Frame of Reference（MRF））来分析，第二部分将采用滑移网格（Sliding Mesh（SM））来作分析。

如上图所示，流过某垂直轴风力机（Vertical Axis Wind Turbine（VAWT））的均匀来流速度为V=10m/s，VAWT直径为12cm且有3个均布叶片，每个叶片弦长为2cm。为简便计，假设其以40RPM的恒定角速度做旋转运动，每个叶片的中心到轮毂中心的距离为0.04m。

实际上，这是对上图所示的三维实际风机的简化二维计算模型。

注意，这是达里厄型垂直轴风力发电机（Darrieus VAWT），其属于升力型风力机。与其对应的是萨沃纽斯型垂直轴风力发电机（Savonius VAWT），这是阻力型风力机。感兴趣的读者可以参考风力机的相关资料。

注意 这里咱们假设风力机的转动是和流体无关的，这和实际情况并不相符，因为流动实际上推动了风力机的旋转运动，即，对于确定的几何形状和质量（转动惯量）的风力机而言，对应外部来流速度的只有唯一的稳定转速。然而，数值模拟由流动所导致的几何体的运动是非常困难的事情，而且需要用到“6DOF solver”（6自由度求解器），首先未必能算，其次算得很慢，再次算得未必靠谱。所以呢，咱们这儿就简单点，同时给定风力机旋转速度和来流的速度，先算个稳态的问题（part1-MRF），再算个瞬态的问题（part2-SM）。

控制方程为运动参考系下的质量守恒和动量守恒（Navier-Stokes方程）。 质量守恒

在以恒定角速度旋转的运动参考系下的Navier-Stokes方程 其中

FLUENT将求解在运动参考系下的上述方程，这个近似所带来的好处就是咱们不再需要处理运动网格了。

注意咱们要求解湍流问题，那么就用 k − ϵ k-\epsilon k−ϵ模型好了。

创建的流体计算域应该比风力机的几何结构大数倍才好，该计算域正是风力机所扰动（波及，影响）到的流场区域。这里用一个外部大圆来作为整个计算域，实际上随便用个几何形状来做为“远场”区域都可以，但是用圆的话相对比较简单些。 边界条件为：

FLUENT将采用有限体积方法求解，先把控制域划分成多个控制体（或单元），然后在每个单元上去积分控制方程，保证每个单元内质量、动量、能量的生成量和通过边界流入的通量是相互平衡的，从而将原本的偏微分控制方程组转化为控制单元的非线性代数方程组，再将其转化为线性化方程组，接下来用迭代方法去求解这些方程组，直到残差满足收敛指标为止。

速度、压力、角速度和湍动能k将参与计算，而其他的值可以通过后处理由这些值推出来，比如剪切力、涡量场等。

CFD的理论相当复杂，具体算法可参考计算流体动力学的相关书籍。

TSR由下式定义： λ = 叶 尖 处 的 速 度 来 流 风 速 = r × ω V \lambda=\frac{叶尖处的速度}{来流风速}=\frac{r\times\omega}{V} λ=来流风速叶尖处的速度​=Vr×ω​ 其中 r = 0.04 m r=0.04m r=0.04m为叶片中心到旋转中心的距离； ω = 40 r p m = 4.1888 r a d / s \omega=40rpm=4.1888rad/s ω=40rpm=4.1888rad/s为旋转角速度； V = 10 m / s V=10m/s V=10m/s为来流速度，这样，叶片中心处的预期速度应该是 r × ω = 0.1676 m / s r\times\omega=0.1676m/s r×ω=0.1676m/s，TSR应该是 0.01676 0.01676 0.01676。

还要引入下面的假设：

新建文件夹vawt_MRF，咱们把后面的算例都放到这个文件夹里面去。

几何模型的构建比较繁琐，这里就不再介绍了，直接把网格给出来。

网格剖分也比较繁琐，这里略去，直接给出网格文件，点击下面网址下载网格文件mesh.msh。

打开FLUENT软件：

选择2D求解器，双精度Double Precision选项，保持默认单核串行就好了，这个算例相对简单，不需要做并行计算。

File -> Read -> Mesh，找到刚才下载的mesh.msh文件，把它读进来。

Setting Up Domain -> Info -> Size，查看下网格单元、面、节点数目。

可以发现这个网格的节点和单元数目并不多，那么计算量也不大了。

Check，检查网格，注意看又没有负体积，计算域的x、y、z上下限是不是跟之前建模的一样，如果不一致，是不是差了一定倍数，那就需要做Scale缩放处理。这里一切都OK的，表明网格没问题。

接下来创建不同区域间的交界面，为了研究网格的旋转，这里定义了多个区域，有叶片区域、轮毂区域和远场区域。

网格划分时用的是“不相符”方法，即，在不同区域的交界面上网格节点不是相互对应吻合的。所以需要告诉FLUENT穿过这些交界面的邻近单元享用同样的信息。

点亮“Boundary Conditions”，点亮“blade_bot_in”，将其Type从“Wall”改变为“Interface”，OK仍旧保持其原有命名。 你会发现左边Tree中多了一项“Mesh Interfaces”，后面会用它来匹配不同区域的交界面。

重复该操作，把如下的边界面全设置成interface类型： blade_bot_in，blade_bot_outer；blade_right_in，blade_right_outer； blade_top_in，blade_top_outer；hub_inner，hub_outer

接下来把这些界面作配对，点亮“Mesh Interfaces”，点击Create/Edit…，

在Create/Edit Mesh Interfaces对话框中，于Mesh Interface选项下写入交界面的名字“blade_bot”，然后在Interface Zone Side 1下选择blade_bot_in，在Interface Zone Side 2下选择blade_bot_outer，最后点击Create，即创建好了blade_bot_in和blade_bot_outer匹配的交界面blade_bot。 重复该操作创建所有4组交界面的配对

这时候，创建好的四个交界面配对就显示出来了

Setup -> General，保持默认Pressure-Based Solver压力基求解器（适用于不可压缩流动）以及Steady稳态问题。

点开Units…设置单位，将angular-velocity单位修改为rpm。

Setup -> Models -> Viscous - Laminar，Edit…， Model选择k-epsilon (2 eqn)，将k-epsilon Model选为Realizable，Ok。 注1：k-epsilon模型中，从理论上来讲Standard、RNG、Realizable是越来越好的。 注2：壁面函数处理方法上，由于这里网格比较粗糙，所以直接用默认的标准壁面函数就好，不要用Enhanced Wall Treatment。

Setup -> Materials -> Fluid -> Air ，保持默认的空气物性参数即可。 密度 ρ = 1.225 k g / m 3 \rho=1.225kg/m^3 ρ=1.225kg/m3，动力粘性系数 μ = 1.7894 × 1 0 − 5 k g / ( m ⋅ s ) \mu=1.7894\times10^{-5}kg/(m\cdot s) μ=1.7894×10−5kg/(m⋅s)。

点亮“Cell Zone Conditions”，注意有5个区域：

点亮fluid-surface_body，点击Edit…，

这个区域的网格是静止的，所以不需要做任何设置，只需要确保Material Name是默认的air就好，OK关掉。

点亮inner，点击Edit…，

勾选Frame Motion，来设置该网格区域的运动参考系，接下来设置角速度和旋转中心。

建模的时候，把风力机的旋转中心放在了坐标原点，所以这里只要确保Rotation-Axis Origin是默认的(0, 0)点就好。

风力机是以40rpm的速度旋转的，将Rotational Velocity设置为40rpm，即可。注意，Relative Specification中是默认的Relative to Cell Zone为absolute，所以这里的转速是绝对速度。

OK确认设置。

点亮blade_top，点击Edit…，

再次勾选Frame Motion，

为了保证叶片区域的旋转和轮毂区域的旋转是一致的，将Relative Specification的Relative to Cell Zone选择为inner，两者相对速度为0。

指定旋转中心的时候要比较小心，叶片中心位于距离风力机旋转轴0.04m的地方，而上方叶片位于120度方位角处（从x轴开始算起），所以旋转中心的位置是：

x = 0.04 × c o s ( 120 ° ) = − 0.02 m x=0.04\times cos(120\degree)=-0.02m x=0.04×cos(120°)=−0.02m y = 0.04 × s i n ( 120 ° ) = 0.034641 m y=0.04\times sin(120\degree)=0.034641m y=0.04×sin(120°)=0.034641m

把该坐标值输入Rotation-Axis Origin (Relative)，OK即可。 同样的操作，将right和bottom叶片网格域也设置好运动方式，这三个叶片网格域的旋转中心汇总如下：

Setup -> Boundary Conditions，点亮farfield1，将其类型从wall变为velocity-inlet，

在Velocity Inlet设置窗口，将Velocity Specification Method改为Components，设置X-Velocity (m/s)为10，

在Turbulence栏，将Specification Method设为Intensity and Length Scale，设置Turbulence Intensity (%)为5，将Turbulent Length Scale (m)设置为1。（这里设置不好的话收敛性很差的……）

OK关闭。

将farfield2类型从wall变为pressure-outlet，

在Pressure Outlet设置窗口，保持Gauge Pressure (pascal)为默认的0，添加和速度进口一样的湍流条件。

Setup -> Boundary Conditions，点亮wall_blade_bot，确认其类型是wall，不然将其改为Wall。

Edit…打开设置窗口，将Wall Motion由Stationary Wall改为Moving Wall，将Motion类型修改为Rotational，并将Rotation-Axis Origin给成和前面blade_bot区域一样的值，即x=-0.02，y=-0.034641，OK。

同样来设置wall_blade_top和wall_blade_right，注意旋转中心坐标是不同的。

可以保存下cas，接下来咱们设置求解器。

Solution -> Solution Methods，设置Scheme为Coupled，其余保持默认。

Solution -> Monitors > Residuals - Print, Plot，Edit…， 将连续方程、x和y速度、k、epsilon的收敛指标全部设置为1e-6，OK。

Create > Moment…，设置监测力矩项。 在Moment Monitor窗口，勾选Options下的Print to Console和Plot，点亮Wall Zones中的wall_blade_right，将Moment Center改为X=0.04m和Y=0m，OK。 同样的操作，监测wall_blade_top上的力矩，注意其Moment Center为X=-0.02m和Y=0.034641m。

Solution -> Solution Initialization， 选择Hybrid Initialization，然后点击initialize来完成初始化流场。 Solution -> Run Calculation，将Number of Iterations设为600，然后点击Calculate开始计算。

计算还是蛮快的，大约465步就收敛了，如果没收敛的话，那么再多算几步计算到收敛吧……

可以看到，收敛曲线还是蛮不错的，计算出的力矩值也区域平稳，表明计算到了稳定收敛状态。

Results -> Graphics -> Contours， Set up…，

勾选Options下的Filled，选择Contours of下的Velocity和Velocity Magnitude，将Levels设置为100，Display。

注意当叶片垂直于流动方向时（最右边那个叶片），会有一个巨大的回流区域，这会给置于该风力机下游的其它风力机造成极大的不利影响，这在设计排布多个VAWT的时候是需要特别注意的。

这里实际上只是旋转风力机的某个瞬时状态，非定常分析将会在后续教程中详解，还会给出非定常的动画展示。

通过对风力机下游远场边界的查看可以发现，流场依旧没有恢复到其自由来流的状态，黄色区域大约是9m/s的样子。

Results -> Graphics -> Contours， Set up…，

勾选Options下的Filled，选择Contours of下的Velocity和Vorticity Magnitude，将Levels设置为100，取消勾选Auto Range，将Min(1/s)和Max(1/s)分别设置为0和2000，取消勾选Clip to Range，Display。

可见涡从叶片尖端释放出来，回想一下一个垂直于流动方向的平板，其流场也是呈现这种特性，即涡从平板的尖端生成。

对于和流动方向近乎垂直的叶片，涡并没有发展得太过剧烈（强度较弱）。

也可以画出很酷的压力云图（它跟速度云图是有关连的哦）和湍动能图（它跟涡量云图是有关系的呢）来，注意把湍动能的尺度设置为Min=0，Max=7J/kg那么出来的图很漂亮的，留给感兴趣的读者自己做吧。

Results -> Reports -> Fluxes，Set Up…

保持Options为默认的Mass Flow Rate，选择Boundaries中的farfield1和farfield2，Compute。

从进口farfield1流入的质量流量为73.5kg/s，从出口farfield2流出的质量流量为73.5kg/s（流入-73.5kg/s），整个计算域的不平衡质量流量为1.65e-8kg/s，这是由计算误差引起的，相对误差大约为2e-10，几乎为0，那么整个计算域的质量流量是守恒的。

同样的方法，可以查看轮毂区域，即风力机区域的质量流量情况，这次选择的是hub_inner和hub_outer，注意这俩实际上是同一个环，只是分属于不同区域的边界，所以两个的值应该是一样的（符号不同是因为它们的法方向不一样）。由于是个包围风力机的封闭圆环边界，所以其质量流量应该是0，那么算给出的值是1e-9，几乎是0了。

Results -> Plots -> XY Plot，Set Up…，

Y Axis Function选为Velocity和Y Velocity，Surfaces点亮wall_blade_right，Plot。 得到了wall_blade_right上y方向的速度，对于最右边的叶片而言，其中点x=0.04m处的y方向速度Vy即其绝对速度（因为右侧叶片位于水平位置，其中心处的旋转速度（牵连速度）是沿着y方向的，其相对速度也是沿着y方向的，所以其绝对速度也是沿着y方向的），那么从图中可以看到x=0.04m处的Vy位于0.165和0.170m/s之间，咱们之前算过，这个点处的绝对速度为0.1676m/s，两者是一致的，表明咱们前面对于运动网格域的施加方法是正确的。

计算TSR的话，是要除上来流速度10m/s的，这样TSR=0.01676。

看下力矩系数（Cm）你会发现它不是0，这意味着有额外的力矩施加在风力机上。没错，你甚至可以用它来计算功率系数（Cp）。

能不能找到一个稳态角速度，当风力机在这个速度下旋转时，不产生额外的力矩呢？即，如果它转快了，那么会有负向力矩使它慢下来；如果它转慢了，有正向力矩让它加速旋转。

为了使用运动参考系方法来找到这样一个速度，可以通过改变给定的角速度值来实现，即Cell Zone Condition中，打开inner并改变它的角速度值直到Cm是0为止。

你会发现稳态角速度大约是80RPM。

当然可以通过加密网格来提高精度，点击下载加密网格文件mesh.msh。